This report fits bivariate Bayesian model similar to those described in Riley et al (2008), Reitsma et al (2005) and Owen et al (2018) in order to get a meta-analysis estimate of the operational characteristics of diagnostic tests for liver fibrosis and cirrhosis based on different biomarkers.
Below we define 4 models, corresponding to those from Riley et al (2008). Each of those models would be fit several times: for each combination of reference diagnostic and experimental diagnostic.
datFile<-"../2021-06-24_RevisedData/HEPSANET_biomarker_cleaned_data.csv"
dat<-read.csv(datFile)
dat<-dat %>%
mutate(ReasontotestHBVSCLF=gsub(pattern=" ",replacement="",case_when(ReasontotestHBVSCLF==""~"N",ReasontotestHBVSCLF=="I"~"L",TRUE~ReasontotestHBVSCLF))) %>%
mutate(
testReason=factor(ReasontotestHBVSCLF,levels=c("C","S","L","F","N")),
studyType=factor(hosp_com)
)
setOfDummyVars<-model.matrix(~studyType+testReason,data=dat)
dat<-cbind(dat,setOfDummyVars[,-1])
The file ../2021-06-24_RevisedData/HEPSANET_biomarker_cleaned_data.csv contains 2990 rows / observations and 91 columns / variables.
Notes (as per Alex’s email from 24 June 2021):
2 samples (SN030, SN021) has reason for testing listed as I, These have been replaced these with L.
Only 2 samples (H014, H009) were the reason for testing is F. These have been merged with S?
dat$testReason[dat$Nomerged %in% c("H014","H009")]<-"S"
dat<-dat %>%
mutate(testReason=factor(as.character(testReason),levels=c("C","S","L","N")))
We filter out participants who are positive for hepatitis C and/or D infection using the variables antihcv, hcvrna, and antihdv (all participants positive for any of these are filtered out). Participants with no records for these variables are retained.
nBefore<-nrow(dat)
dat<-dat %>%
filter((is.na(antihcv) | antihcv!=1) &
(is.na(hcvrna) | hcvrna!=1) &
(is.na(antihdv) |antihdv!=1))
This has filtered out 45 participants, leaving rnrow(dat)` participants for the analysis.
sumFun<-function(dat,vars,type,digs=2,na.rm=TRUE){
# dat = data frame to be summarised
# vars = character vector with variables names
# type = character vector of same length as vars specifying the type of summary to compute (needs to be one of "mean_sd","median_iqr","n_perc")
# digs = integer specifying the number of decimal digits to report (defaults to 2)
# na.rm = whether to remove NAs or not (defaults to TRUE)
res<-data.frame(
var=character(0),
stat1=character(0),
stat2=character(0)
)
for(i in 1:length(vars)){
if(type[i]=="mean_sd"){
stat1<-format(nsmall=digs,round(digits=digs,mean(dat[,vars[i]],na.rm=T)))
stat2<-paste(sep="","(",format(nsmall=digs,round(digits=digs,sd(dat[,vars[i]],na.rm=T))),")")
#varName<-paste(sep=" ",vars[i],"; mean (sd)")
varName<-paste(sep=" ","mean (sd)")
res<-rbind(res,c(varName,stat1,stat2))
rm(varName,stat1,stat2)
}else if(type[i]=="median_iqr"){
stat1<-format(nsmall=digs,round(digits=digs,median(dat[,vars[i]],na.rm=T)))
stat2<-paste(sep="","(",paste(collapse=",",format(nsmall=digs,round(digits=digs,quantile(dat[,vars[i]],probs=c(0.25,0.75),na.rm=T)))),")")
#varName<-paste(sep=" ",vars[i],"; median (IQR)")
varName<-paste(sep=" ","median (IQR)")
res<-rbind(res,c(varName,stat1,stat2))
rm(varName,stat1,stat2)
}else if(type[i]=="n_perc"){
for(val in sort(unique(dat[,vars[i]]))){
stat1<-sum(dat[,vars[i]]==val,na.rm=T)
stat2<-paste(sep="","(",format(nsmall=digs,round(digits=digs,100*sum(dat[,vars[i]]==val,na.rm=T)/sum(!is.na(dat[,vars[i]])))),"%)")
#varName<-paste(sep="",vars[i],"; ",val,"; n (%)")
varName<-paste(sep="",val,"; n (%)")
res<-rbind(res,c(varName,stat1,stat2))
rm(varName,stat1,stat2)
}
}
}
colnames(res)<-c("variable","stat1","stat2")
return(res)
}
datSum<-sumFun(dat,
vars=c("hosp_com","site","sex","age","bmi","te","apri","gpr","fib4","alt"),
type=c("n_perc","n_perc","n_perc","median_iqr","mean_sd","median_iqr","median_iqr","median_iqr","median_iqr","median_iqr","median_iqr"))
datSum %>%
kable(caption="Patient characteristics",row.names=F) %>%
kable_styling(full_width=F) %>%
pack_rows(start_row=1,end_row=2,group_label=paste(sep="","Type of study; n = ",sum(!is.na(dat$hosp_com))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$hosp_com))/nrow(dat))),"%)")) %>%
pack_rows(start_row=3,end_row=13,group_label=paste(sep="","Site; n = ",sum(!is.na(dat$site))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$site))/nrow(dat))),"%)")) %>%
pack_rows(start_row=14,end_row=15,group_label=paste(sep="","Sex; n = ",sum(!is.na(dat$sex))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$age))/nrow(dat))),"%)")) %>%
pack_rows(start_row=16,end_row=16,group_label=paste(sep="","Age; n = ",sum(!is.na(dat$age))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$age))/nrow(dat))),"%)")) %>%
pack_rows(start_row=17,end_row=17,group_label=paste(sep="","BMI; n = ",sum(!is.na(dat$bmi))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$bmi))/nrow(dat))),"%)")) %>%
pack_rows(start_row=18,end_row=18,group_label=paste(sep="","Transient elastography; n = ",sum(!is.na(dat$te))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$te))/nrow(dat))),"%)")) %>%
pack_rows(start_row=19,end_row=19,group_label=paste(sep="","APRI (AST to platelet ratio); n = ",sum(!is.na(dat$apri))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$apri))/nrow(dat))),"%)")) %>%
pack_rows(start_row=20,end_row=20,group_label=paste(sep="","GPR (GGT to platelet ratio); n = ",sum(!is.na(dat$gpr))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$gpr))/nrow(dat))),"%)")) %>%
pack_rows(start_row=21,end_row=21,group_label=paste(sep="","Fibrosis 4 score; n = ",sum(!is.na(dat$bmi))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$bmi))/nrow(dat))),"%)")) %>%
pack_rows(start_row=22,end_row=22,group_label=paste(sep="","ALT; n = ",sum(!is.na(dat$bmi))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$bmi))/nrow(dat))),"%)"))
| variable | stat1 | stat2 |
|---|---|---|
| Type of study; n = 2945 (100.0%) | ||
| Community; n (%) | 735 | (24.96%) |
| Hospital; n (%) | 2210 | (75.04%) |
| Site; n = 2945 (100.0%) | ||
| Burkina Faso; n (%) | 35 | (1.19%) |
| Cape Town; n (%) | 155 | (5.26%) |
| Ethiopia; n (%) | 1047 | (35.55%) |
| Gambia; n (%) | 782 | (26.55%) |
| Malawi; n (%) | 100 | (3.40%) |
| Nigeria; n (%) | 190 | (6.45%) |
| Senegal_HPB; n (%) | 84 | (2.85%) |
| Senegal_Maud; n (%) | 275 | (9.34%) |
| Senegal_PRF; n (%) | 151 | (5.13%) |
| Stellenbosch; n (%) | 36 | (1.22%) |
| Zambia; n (%) | 90 | (3.06%) |
| Sex; n = 2945 (100.0%) | ||
| Female; n (%) | 1138 | (38.64%) |
| Male; n (%) | 1807 | (61.36%) |
| Age; n = 2945 (100.0%) | ||
| median (IQR) | 34.00 | (28.00,41.00) |
| BMI; n = 2583 (87.7%) | ||
| mean (sd) | 23.12 | (4.71) |
| Transient elastography; n = 2945 (100.0%) | ||
| median (IQR) | 5.60 | (4.50,7.10) |
| APRI (AST to platelet ratio); n = 2945 (100.0%) | ||
| median (IQR) | 0.32 | (0.22,0.48) |
| GPR (GGT to platelet ratio); n = 2261 (76.8%) | ||
| median (IQR) | 0.12 | (0.07,0.18) |
| Fibrosis 4 score; n = 2583 (87.7%) | ||
| median (IQR) | 0.85 | (0.58,1.30) |
| ALT; n = 2583 (87.7%) | ||
| median (IQR) | 25.00 | (19.00,35.00) |
Reference standards (which are assumed to be gold standards) are based on variable te, transient elastrography result. Specifically:
In the dataset, these 3 gold standards are encoded by variables te_115, te_95 and te_79.
dat<-dat %>%
mutate(
te_115=ifelse(te>11.5,1,0),
te_95=ifelse(te>9.5,1,0),
te_79=ifelse(te>7.9,1,0),
)
apri - AST to platelet ratio index. For F4 use threshold of \(\geq 2.0\) and for F2 use \(\geq 1.5\). In the dataset, these 2 experimental diagnostics are encoded by variables apri_20 and apri_15. Note: we will also derive an optimum threshold.
gpr - GGT to platelet ratio. For F4 use threshold \(\geq 0.56\) and for F2 use \(\geq 0.32\). In the dataset, these 2 experimental diagnostics are encoded by variables gpr_56 and gpr_32.
fib4 - fibrosis 4 score. We need to derive a score for this; no a prior defined threshold.
alt - ALT. We need to derive a score for this; no a prior defined threshold.
dat<-dat %>%
mutate(
apri_20=ifelse(apri>=2,1,0),
apri_15=ifelse(apri>=1.5,1,0),
gpr_56=ifelse(gpr>=0.56,1,0),
gpr_32=ifelse(gpr>=0.32,1,0),
)
Let’s just assume we could analyse all individual-level data as is, simply pulling out the sensitivity, specificity rather than fitting a model that considers random effects and individual- and study-level covariates.
gr<-rbind(
c(11.5,"2.0","0.56"),
c(9.5,"2.0","0.56"),
c(7.9,"1.5","0.32")
)
rnames<-character(0)
for(s in c("All",sort(unique(dat$Sitecountry)))){
for(i in 1:3){
rnames<-c(rnames,paste(sep="_",paste(sep="","Site",s),paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))))
}
}
rnames<-gsub(rnames,pattern="SiteAll",replacement="All")
sensSpecMat<-matrix(nrow=3*(length(unique(dat$Sitecountry))+1),ncol=8)
rownames(sensSpecMat)<-rnames
colnames(sensSpecMat)<-c("site","te","apri","gpr","apri_sens","apri_spec","gpr_sens","gpr_spec")
sensSpecMat<-as.data.frame(sensSpecMat)
for(s in c("All",sort(unique(dat$Sitecountry)))){
if(s=="All"){
datTmp<-dat
}else{
datTmp<-dat %>% filter(Sitecountry==s)
s<-paste(sep="","Site",s)
}
for(i in 1:nrow(gr)){
gthr<-gr[i,1]
thr1<-gr[i,2]
thr2<-gr[i,3]
sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"site"]<-gsub(s,pattern="Site",replacement="")
sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"te"]<-gthr
sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"apri"]<-thr1
sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"gpr"]<-thr2
sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"apri_sens"]<-mean(datTmp[datTmp[,paste(sep="_","te",gsub(gthr,pattern="\\.",replacement=""))]==1,paste(sep="_","apri",gsub(thr1,pattern="\\.",replacement=""))])
sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"gpr_sens"]<-mean(datTmp[datTmp[,paste(sep="_","te",gsub(gthr,pattern="\\.",replacement=""))]==1,paste(sep="_","gpr",gsub(thr2,pattern="0\\.",replacement=""))],na.rm=T)
sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"apri_spec"]<-mean(1-datTmp[datTmp[,paste(sep="_","te",gsub(gthr,pattern="\\.",replacement=""))]==0,paste(sep="_","apri",gsub(thr1,pattern="\\.",replacement=""))])
sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"gpr_spec"]<-mean(1-datTmp[datTmp[,paste(sep="_","te",gsub(gthr,pattern="\\.",replacement=""))]==0,paste(sep="_","gpr",gsub(thr2,pattern="0\\.",replacement=""))],na.rm=T)
}
}
sensSpecMatKable<-sensSpecMat
for(i in 1:nrow(sensSpecMat)){
for(j in 5:ncol(sensSpecMat)){
if(!is.na(sensSpecMat[i,j])){
sensSpecMatKable[i,j]<-paste(sep="",format(nsmall=1,round(digits=1,100*sensSpecMat[i,j])),"%")
}else{
sensSpecMatKable[i,j]<-"-"
}
}
}
rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te11.5_apri2.0_gpr0.56")]<-gsub(rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te11.5_apri2.0_gpr0.56")],pattern="_te11.5_apri2.0_gpr0.56",replacement="; te > 11.5 (apri > 2.0, gpr > 0.56)")
rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te9.5_apri2.0_gpr0.56")]<-gsub(rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te9.5_apri2.0_gpr0.56")],pattern="_te9.5_apri2.0_gpr0.56",replacement="; te > 9.5 (apri > 2.0, gpr > 0.56)")
rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te7.9_apri1.5_gpr0.32")]<-gsub(rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te7.9_apri1.5_gpr0.32")],pattern="_te7.9_apri1.5_gpr0.32",replacement="; te > 7.9 (apri > 1.5, gpr > 0.32)")
sensSpecMatKable <- cbind("",sensSpecMatKable)
sensSpecMatKable %>%
dplyr::select(!site) %>%
kable(caption="Naive sensitivity and specificity estimates (expressed as percentages).",col.names=c(" ","TE","APRI","GPR","Sensitivity","Specificity","Sensitivity","Specificity"),row.names = FALSE) %>%
kable_styling(full_width = FALSE) %>%
add_header_above(c(" "=1,"Diagnostic thresholds"=3,"APRI"=2,"GPR"=2)) %>%
pack_rows(start_row=1,end_row=3,group_label="All samples") %>%
pack_rows(start_row=4,end_row=6,group_label="Site 1") %>%
pack_rows(start_row=7,end_row=9,group_label="Site 2") %>%
pack_rows(start_row=10,end_row=12,group_label="Site 3") %>%
pack_rows(start_row=13,end_row=15,group_label="Site 4") %>%
pack_rows(start_row=16,end_row=18,group_label="Site 5") %>%
pack_rows(start_row=19,end_row=21,group_label="Site 6") %>%
pack_rows(start_row=22,end_row=24,group_label="Site 7") %>%
pack_rows(start_row=25,end_row=27,group_label="Site 8") %>%
pack_rows(start_row=28,end_row=30,group_label="Site 9") %>%
pack_rows(start_row=31,end_row=33,group_label="Site 10") %>%
pack_rows(start_row=34,end_row=36,group_label="Site 11")
| TE | APRI | GPR | Sensitivity | Specificity | Sensitivity | Specificity | |
|---|---|---|---|---|---|---|---|
| All samples | |||||||
| 11.5 | 2.0 | 0.56 | 20.2% | 99.5% | 31.6% | 98.0% | |
| 9.5 | 2.0 | 0.56 | 15.5% | 99.5% | 26.5% | 98.3% | |
| 7.9 | 1.5 | 0.32 | 14.4% | 99.2% | 37.1% | 94.6% | |
| Site 1 | |||||||
| 11.5 | 2.0 | 0.56 | 10.5% | 99.8% | 18.6% | 99.7% | |
| 9.5 | 2.0 | 0.56 | 8.5% | 99.8% | 15.1% | 99.7% | |
| 7.9 | 1.5 | 0.32 | 10.1% | 99.3% | 29.3% | 97.5% | |
| Site 2 | |||||||
| 11.5 | 2.0 | 0.56 | 0.0% | 99.3% | 22.2% | 95.0% | |
| 9.5 | 2.0 | 0.56 | 0.0% | 99.3% | 20.0% | 96.2% | |
| 7.9 | 1.5 | 0.32 | 5.3% | 99.1% | 44.7% | 92.9% | |
| Site 3 | |||||||
| 11.5 | 2.0 | 0.56 |
|
100.0% |
|
100.0% | |
| 9.5 | 2.0 | 0.56 | 0.0% | 100.0% | 0.0% | 100.0% | |
| 7.9 | 1.5 | 0.32 | 0.0% | 100.0% | 20.0% | 98.2% | |
| Site 4 | |||||||
| 11.5 | 2.0 | 0.56 | 70.0% | 100.0% | 83.3% | 98.6% | |
| 9.5 | 2.0 | 0.56 | 65.6% | 100.0% | 78.1% | 98.5% | |
| 7.9 | 1.5 | 0.32 | 63.2% | 100.0% | 84.2% | 95.2% | |
| Site 5 | |||||||
| 11.5 | 2.0 | 0.56 | 26.1% | 100.0% |
|
|
|
| 9.5 | 2.0 | 0.56 | 18.2% | 100.0% |
|
|
|
| 7.9 | 1.5 | 0.32 | 18.9% | 99.3% |
|
|
|
| Site 6 | |||||||
| 11.5 | 2.0 | 0.56 | 40.0% | 100.0% |
|
|
|
| 9.5 | 2.0 | 0.56 | 28.6% | 100.0% |
|
|
|
| 7.9 | 1.5 | 0.32 | 23.1% | 100.0% |
|
|
|
| Site 7 | |||||||
| 11.5 | 2.0 | 0.56 | 33.3% | 97.5% | 66.7% | 86.1% | |
| 9.5 | 2.0 | 0.56 | 26.7% | 98.7% | 66.7% | 90.9% | |
| 7.9 | 1.5 | 0.32 | 25.0% | 98.5% | 44.4% | 76.7% | |
| Site 8 | |||||||
| 11.5 | 2.0 | 0.56 | 21.4% | 99.1% | 42.9% | 97.4% | |
| 9.5 | 2.0 | 0.56 | 10.0% | 99.1% | 27.6% | 97.6% | |
| 7.9 | 1.5 | 0.32 | 4.1% | 98.6% | 31.4% | 91.6% | |
| Site 9 | |||||||
| 11.5 | 2.0 | 0.56 | 13.3% | 100.0% | 33.3% | 98.5% | |
| 9.5 | 2.0 | 0.56 | 11.8% | 100.0% | 29.4% | 98.5% | |
| 7.9 | 1.5 | 0.32 | 10.0% | 100.0% | 50.0% | 95.4% | |
| Site 10 | |||||||
| 11.5 | 2.0 | 0.56 | 100.0% | 100.0% | 100.0% | 100.0% | |
| 9.5 | 2.0 | 0.56 | 50.0% | 100.0% | 50.0% | 100.0% | |
| 7.9 | 1.5 | 0.32 | 50.0% | 100.0% | 50.0% | 100.0% | |
| Site 11 | |||||||
| 11.5 | 2.0 | 0.56 | 8.7% | 99.2% | 9.1% | 97.5% | |
| 9.5 | 2.0 | 0.56 | 9.4% | 99.6% | 13.3% | 98.3% | |
| 7.9 | 1.5 | 0.32 | 10.2% | 99.5% | 25.0% | 97.1% | |
Question:
Sensitivity and specificity are extremely variable from study to study. I find it difficult to understand this. You would expect some heterogeneity due to how certain diagnostics are delivered, some underlying characteristics that impct test results which vary from population to population, but not quite to this extent. Positive and negative predictive values can vary a lot from one population to another, but that is because these, unlike sensitvity and specificity, depend on the prevalence of the condition under investigation.
Would be good to discuss this.
Earlier versions of this document (20210627 and earlier) compared several different models. From this, the decision was made for a model with individual- and study-level covariates and study-specific random effects
This is model (4) from Riley et al (2008).
This model estimates summary sensitivity and specificity with study random effects and individual- and study-level covariates.
We include the following co-variates:
As noted earlier in this document, there exist recommended thresholds for APRI (>2.0 for F4 and >1.5 for F2) and GPR (>0.56 for F4 and >0.32 for F2). Together with the 2 thresholds for the reference standard for F4 (te>11.5 and te>9.5) and 1 threshold for F2 (te>7.9), this gives 6 models for pre-defined models. The overall sensitivities and specificities for these models ar reported in Tables 3 and 4 below. Model-specific tables giving ORs for the different co-variates are given in the appendix.
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_115==1,1,2)
datJags$y<-dat$apri_20
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4_te115_apri20 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4_te115_apri20<-coda.samples(model=jagsModel4_te115_apri20,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
tmp<-try(summary(parsModel4_te115_apri20)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te115_apri20))
mod_te115_apri20<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_95==1,1,2)
datJags$y<-dat$apri_20
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4_te95_apri20 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4_te95_apri20<-coda.samples(model=jagsModel4_te95_apri20,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
tmp<-try(summary(parsModel4_te95_apri20)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te95_apri20))
mod_te95_apri20<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_79==1,1,2)
datJags$y<-dat$apri_15
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4_te79_apri15 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4_te79_apri15<-coda.samples(model=jagsModel4_te79_apri15,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
tmp<-try(summary(parsModel4_te79_apri15)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te79_apri15))
mod_te79_apri15<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_115==1,1,2)
datJags$y<-dat$gpr_56
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4_te115_gpr56 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4_te115_gpr56<-coda.samples(model=jagsModel4_te115_gpr56,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
tmp<-try(summary(parsModel4_te115_gpr56)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te115_gpr56))
mod_te115_gpr56<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_95==1,1,2)
datJags$y<-dat$gpr_56
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4_te95_gpr56 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 10000) # needs to be large!
parsModel4_te95_gpr56<-coda.samples(model=jagsModel4_te95_gpr56,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
tmp<-try(summary(parsModel4_te95_gpr56)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te95_gpr56))
mod_te95_gpr56<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_79==1,1,2)
datJags$y<-dat$gpr_32
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4_te79_gpr32 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4_te79_gpr32<-coda.samples(model=jagsModel4_te79_gpr32,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
tmp<-try(summary(parsModel4_te79_gpr32)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te79_gpr32))
mod_te79_gpr32<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
save(list=c("jagsModel4_te115_apri20","jagsModel4_te95_apri20","jagsModel4_te79_apri15","jagsModel4_te115_gpr56","jagsModel4_te95_gpr56","jagsModel4_te79_gpr32","parsModel4_te115_apri20","parsModel4_te95_apri20","parsModel4_te79_apri15","parsModel4_te115_gpr56","parsModel4_te95_gpr56","parsModel4_te79_gpr32"),file=paste(sep="","fixedThresholdModelsResultsList_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
## Error in ar.yw.default(x, aic = aic, order.max = order.max, na.action = na.action, :
## NA/NaN/Inf in foreign function call (arg 2)
## Error in ar.yw.default(x, aic = aic, order.max = order.max, na.action = na.action, :
## NA/NaN/Inf in foreign function call (arg 2)
## Error in ar.yw.default(x, aic = aic, order.max = order.max, na.action = na.action, :
## NA/NaN/Inf in foreign function call (arg 2)
## Error in ar.yw.default(x, aic = aic, order.max = order.max, na.action = na.action, :
## NA/NaN/Inf in foreign function call (arg 2)
## Error in ar.yw.default(x, aic = aic, order.max = order.max, na.action = na.action, :
## NA/NaN/Inf in foreign function call (arg 2)
## Error in ar.yw.default(x, aic = aic, order.max = order.max, na.action = na.action, :
## NA/NaN/Inf in foreign function call (arg 2)
## Error in ar.yw.default(x, aic = aic, order.max = order.max, na.action = na.action, :
## NA/NaN/Inf in foreign function call (arg 2)
## Error in ar.yw.default(x, aic = aic, order.max = order.max, na.action = na.action, :
## NA/NaN/Inf in foreign function call (arg 2)
| te>11.5 | te>9.5 | te>7.9 |
|---|---|---|
|
sensitivity (95% CI) |
||
| 28.76% | 14.76% | 7.32% |
| (5.70%,54.83%) | (2.70%,28.52%) | (1.87%,13.93%) |
|
specificity (95% CI) |
||
| 99.18% | 99.16% | 99.16% |
| (98.20%,99.95%) | (98.18%,99.90%) | (98.17%,99.93%) |
| te>11.5 | te>9.5 | te>7.9 |
|---|---|---|
|
sensitivity (95% CI) |
||
| 22.68% | 20.43% | 18.88% |
| (3.33%,47.31%) | (5.34%,37.95%) | (8.12%,30.73%) |
|
specificity (95% CI) |
||
| 98.45% | 98.88% | 96.76% |
| (97.11%,99.56%) | (97.88%,99.69%) | (94.83%,98.50%) |
apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43
apri115ResultsList<-list()
apri115ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(apriVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
apr<-apriVect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_115==1,1,2)
datJags$y<-ifelse(dat$apri>=apr,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(apriVect)){
apr<-apriVect[i]
apri115ModList[[paste(sep="","apri_mod_",apr)]]<-parList[[i]]$model
apri115ModList[[paste(sep="","apri_samples_",apr)]]<-parList[[i]]$pars
apri115ResultsList[[paste(sep="","apri",apr)]]<-parList[[i]]$results
}
rm(parList)
save(apri115ResultsList,apri115ModList,file=paste(sep="","apriResultsList_te115_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43
load("apriResultsList_te115_20210805.RData")
sensVect<-numeric(length(apriVect))
specVect<-numeric(length(apriVect))
for(i in 1:length(apriVect)){
sensVect[i]<-apri115ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-apri115ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
apri=apriVect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-apriVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nAPRI threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for APRI with diagnosis threshold te > 11.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on APRI (from 0 to 13.26 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
apri115ResultsList[[paste(sep="","apri",bestThr)]] %>%
kable(caption="Best (according to Youden's J) APRI model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 0.5804 | (0.0958,1.2771) |
| hazardous alcohol drinking (OR for specificity) | 1.0163 | (0.5602,1.5236) |
| suspected liver disease screening (OR for sensitivity) | 3.7027 | (0.6861,8.2046) |
| suspected liver disease screening (OR for specificity) | 0.1457 | (0.0788,0.2173) |
| hospital-based study (OR for sensitivity) | 2.4776 | (0.2120,6.5763) |
| hospital-based study (OR for specificity) | 1.4508 | (1.0241,1.9267) |
| random effects variance (logit sensitivity) | 1.2545 | (0.1256,3.4243) |
| random effects variance (logit specificity) | 0.7847 | (0.2034,1.6496) |
| overall sensitivity | 0.6222 | (0.3468,0.8748) |
| overall specificity | 0.8564 | (0.8162,0.8962) |
apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43
apri95ResultsList<-list()
apri95ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(apriVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
apr<-apriVect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_95==1,1,2)
datJags$y<-ifelse(dat$apri>=apr,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(apriVect)){
apr<-apriVect[i]
apri95ModList[[paste(sep="","apri_mod_",apr)]]<-parList[[i]]$model
apri95ModList[[paste(sep="","apri_samples_",apr)]]<-parList[[i]]$pars
apri95ResultsList[[paste(sep="","apri",apr)]]<-parList[[i]]$results
}
rm(parList)
save(apri95ResultsList,apri95ModList,file=paste(sep="","apriResultsList_te95_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43
load("apriResultsList_te95_20210805.RData")
sensVect<-numeric(length(apriVect))
specVect<-numeric(length(apriVect))
for(i in 1:length(apriVect)){
sensVect[i]<-apri95ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-apri95ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
apri=apriVect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-apriVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nAPRI threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for APRI with diagnosis threshold te > 9.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on APRI (from 0 to 13.26 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
apri95ResultsList[[paste(sep="","apri",bestThr)]] %>%
kable(caption="Best (according to Youden's J) APRI model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 0.3891 | (0.0839, 0.8028) |
| hazardous alcohol drinking (OR for specificity) | 1.0314 | (0.5904, 1.5065) |
| suspected liver disease screening (OR for sensitivity) | 13.3700 | (3.6185,26.6388) |
| suspected liver disease screening (OR for specificity) | 0.1644 | (0.1054, 0.2328) |
| hospital-based study (OR for sensitivity) | 0.9215 | (0.2037, 1.8660) |
| hospital-based study (OR for specificity) | 1.4098 | (1.0479, 1.7939) |
| random effects variance (logit sensitivity) | 0.7792 | (0.1152, 1.9096) |
| random effects variance (logit specificity) | 0.7465 | (0.1985, 1.5735) |
| overall sensitivity | 0.7049 | (0.5192, 0.8895) |
| overall specificity | 0.7750 | (0.7263, 0.8235) |
apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43
apri79ResultsList<-list()
apri79ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(apriVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
apr<-apriVect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_79==1,1,2)
datJags$y<-ifelse(dat$apri>=apr,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(apriVect)){
apr<-apriVect[i]
apri79ModList[[paste(sep="","apri_mod_",apr)]]<-parList[[i]]$model
apri79ModList[[paste(sep="","apri_samples_",apr)]]<-parList[[i]]$pars
apri79ResultsList[[paste(sep="","apri",apr)]]<-parList[[i]]$results
}
rm(parList)
save(apri79ResultsList,apri79ModList,file=paste(sep="","apriResultsList_te79_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43
load("apriResultsList_te79_20210805.RData")
sensVect<-numeric(length(apriVect))
specVect<-numeric(length(apriVect))
for(i in 1:length(apriVect)){
sensVect[i]<-apri79ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-apri79ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
apri=apriVect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-apriVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nAPRI threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for APRI with diagnosis threshold te > 7.9kPa (F2).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on APRI (from 0 to 13.26 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
apri79ResultsList[[paste(sep="","apri",bestThr)]] %>%
kable(caption="Best (according to Youden's J) APRI model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 0.3985 | (0.1254, 0.7339) |
| hazardous alcohol drinking (OR for specificity) | 0.9596 | (0.5289, 1.4377) |
| suspected liver disease screening (OR for sensitivity) | 10.7007 | (4.5094,18.3832) |
| suspected liver disease screening (OR for specificity) | 0.1582 | (0.0914, 0.2302) |
| hospital-based study (OR for sensitivity) | 0.6865 | (0.2929, 1.1543) |
| hospital-based study (OR for specificity) | 1.3991 | (1.0028, 1.8296) |
| random effects variance (logit sensitivity) | 0.5718 | (0.1087, 1.3040) |
| random effects variance (logit specificity) | 0.9175 | (0.2338, 1.9148) |
| overall sensitivity | 0.6674 | (0.5289, 0.8061) |
| overall specificity | 0.8091 | (0.7597, 0.8563) |
Note that the GPR data is the only one of APRI, GPR, ALT, FIB4 to have missing values. These are ignored when selecting the threshold points. The Bayesian models will integrate them out in a principled fashion.
gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40
gpr115ResultsList<-list()
gpr115ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(gprVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
gpr<-gprVect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_115==1,1,2)
datJags$y<-ifelse(dat$gpr>=gpr,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(gprVect)){
gpr<-gprVect[i]
gpr115ModList[[paste(sep="","gpr_mod_",gpr)]]<-parList[[i]]$model
gpr115ModList[[paste(sep="","gpr_samples_",gpr)]]<-parList[[i]]$pars
gpr115ResultsList[[paste(sep="","gpr",gpr)]]<-parList[[i]]$results
}
rm(parList)
save(gpr115ResultsList,gpr115ModList,file=paste(sep="","gprResultsList_te115_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40
load("gprResultsList_te115_20210805.RData")
sensVect<-numeric(length(gprVect))
specVect<-numeric(length(gprVect))
for(i in 1:length(gprVect)){
sensVect[i]<-gpr115ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-gpr115ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
gpr=gprVect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-gprVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nGPR threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for GPR with diagnosis threshold te > 11.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on GPR (from 0 to 13.65 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
gpr115ResultsList[[paste(sep="","gpr",bestThr)]] %>%
kable(caption="Best (according to Youden's J) GPR model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 1.1523 | (0.0577, 3.2957) |
| hazardous alcohol drinking (OR for specificity) | 0.5383 | (0.3125, 0.8002) |
| suspected liver disease screening (OR for sensitivity) | 11.3094 | (0.3977,34.4099) |
| suspected liver disease screening (OR for specificity) | 0.3301 | (0.2145, 0.4602) |
| hospital-based study (OR for sensitivity) | 1.8280 | (0.0087, 5.5713) |
| hospital-based study (OR for specificity) | 0.4535 | (0.3272, 0.5837) |
| random effects variance (logit sensitivity) | 1.1528 | (0.0905, 3.3541) |
| random effects variance (logit specificity) | 0.4853 | (0.1088, 1.0596) |
| overall sensitivity | 0.8622 | (0.6569, 0.9994) |
| overall specificity | 0.6979 | (0.6299, 0.7637) |
gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40
gpr95ResultsList<-list()
gpr95ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(gprVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
gpr<-gprVect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_95==1,1,2)
datJags$y<-ifelse(dat$gpr>=gpr,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(gprVect)){
gpr<-gprVect[i]
gpr95ModList[[paste(sep="","gpr_mod_",gpr)]]<-parList[[i]]$model
gpr95ModList[[paste(sep="","gpr_samples_",gpr)]]<-parList[[i]]$pars
gpr95ResultsList[[paste(sep="","gpr",gpr)]]<-parList[[i]]$results
}
rm(parList)
save(gpr95ResultsList,gpr95ModList,file=paste(sep="","gprResultsList_te95_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40
load("gprResultsList_te95_20210805.RData")
sensVect<-numeric(length(gprVect))
specVect<-numeric(length(gprVect))
for(i in 1:length(gprVect)){
sensVect[i]<-gpr95ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-gpr95ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
gpr=gprVect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-gprVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nGPR threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for GPR with diagnosis threshold te > 9.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on GPR (from 0 to 13.65 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
gpr95ResultsList[[paste(sep="","gpr",bestThr)]] %>%
kable(caption="Best (according to Youden's J) GPR model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 1.5598 | (0.1637, 3.8903) |
| hazardous alcohol drinking (OR for specificity) | 0.5464 | (0.3089, 0.8106) |
| suspected liver disease screening (OR for sensitivity) | 9.6354 | (2.1353,20.9855) |
| suspected liver disease screening (OR for specificity) | 0.3326 | (0.2094, 0.4632) |
| hospital-based study (OR for sensitivity) | 0.9495 | (0.0844, 2.3562) |
| hospital-based study (OR for specificity) | 0.4307 | (0.3103, 0.5705) |
| random effects variance (logit sensitivity) | 0.8032 | (0.0865, 2.0909) |
| random effects variance (logit specificity) | 0.5040 | (0.1185, 1.0852) |
| overall sensitivity | 0.8295 | (0.6568, 0.9739) |
| overall specificity | 0.7142 | (0.6402, 0.7813) |
gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40
gpr79ResultsList<-list()
gpr79ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(gprVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
gpr<-gprVect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_79==1,1,2)
datJags$y<-ifelse(dat$gpr>=gpr,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(gprVect)){
gpr<-gprVect[i]
gpr79ModList[[paste(sep="","gpr_mod_",gpr)]]<-parList[[i]]$model
gpr79ModList[[paste(sep="","gpr_samples_",gpr)]]<-parList[[i]]$pars
gpr79ResultsList[[paste(sep="","gpr",gpr)]]<-parList[[i]]$results
}
rm(parList)
save(gpr79ResultsList,gpr79ModList,file=paste(sep="","gprResultsList_te79_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40
load("gprResultsList_te79_20210805.RData")
sensVect<-numeric(length(gprVect))
specVect<-numeric(length(gprVect))
for(i in 1:length(gprVect)){
sensVect[i]<-gpr79ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-gpr79ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
gpr=gprVect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-gprVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nGPR threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for GPR with diagnosis threshold te > 7.9kPa (F2).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on GPR (from 0 to 13.65 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
gpr79ResultsList[[paste(sep="","gpr",bestThr)]] %>%
kable(caption="Best (according to Youden's J) GPR model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 0.7777 | (0.1482, 1.6472) |
| hazardous alcohol drinking (OR for specificity) | 0.5182 | (0.2656, 0.7933) |
| suspected liver disease screening (OR for sensitivity) | 7.0267 | (2.6486,12.4961) |
| suspected liver disease screening (OR for specificity) | 0.3135 | (0.1959, 0.4376) |
| hospital-based study (OR for sensitivity) | 1.7557 | (0.3101, 3.6954) |
| hospital-based study (OR for specificity) | 0.4872 | (0.3401, 0.6424) |
| random effects variance (logit sensitivity) | 0.8070 | (0.1077, 2.0433) |
| random effects variance (logit specificity) | 0.5383 | (0.1208, 1.1604) |
| overall sensitivity | 0.8364 | (0.7037, 0.9576) |
| overall specificity | 0.6119 | (0.5305, 0.6913) |
altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40
alt115ResultsList<-list()
alt115ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(altVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
alt<-altVect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_115==1,1,2)
datJags$y<-ifelse(dat$alt>=alt,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(altVect)){
alt<-altVect[i]
alt115ModList[[paste(sep="","alt_mod_",alt)]]<-parList[[i]]$model
alt115ModList[[paste(sep="","alt_samples_",alt)]]<-parList[[i]]$pars
alt115ResultsList[[paste(sep="","alt",alt)]]<-parList[[i]]$results
}
rm(parList)
save(alt115ResultsList,alt115ModList,file=paste(sep="","altResultsList_te115_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40
load("altResultsList_te115_20210805.RData")
sensVect<-numeric(length(altVect))
specVect<-numeric(length(altVect))
for(i in 1:length(altVect)){
sensVect[i]<-alt115ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-alt115ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
alt=altVect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-altVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nALT threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for ALT with diagnosis threshold te > 11.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on ALT (from 1 to 198 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
alt115ResultsList[[paste(sep="","alt",bestThr)]] %>%
kable(caption="Best (according to Youden's J) ALT model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 0.9139 | (0.1596,1.9809) |
| hazardous alcohol drinking (OR for specificity) | 0.7295 | (0.4582,1.0217) |
| suspected liver disease screening (OR for sensitivity) | 1.0797 | (0.2004,2.2190) |
| suspected liver disease screening (OR for specificity) | 0.2075 | (0.1458,0.2730) |
| hospital-based study (OR for sensitivity) | 1.9812 | (0.1984,4.9471) |
| hospital-based study (OR for specificity) | 0.6845 | (0.5068,0.8671) |
| random effects variance (logit sensitivity) | 0.4747 | (0.0666,1.1553) |
| random effects variance (logit specificity) | 0.8942 | (0.2354,1.8993) |
| overall sensitivity | 0.7080 | (0.4767,0.9099) |
| overall specificity | 0.7311 | (0.6768,0.7824) |
altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40
alt95ResultsList<-list()
alt95ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(altVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
alt<-altVect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_95==1,1,2)
datJags$y<-ifelse(dat$alt>=alt,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(altVect)){
alt<-altVect[i]
alt95ModList[[paste(sep="","alt_mod_",alt)]]<-parList[[i]]$model
alt95ModList[[paste(sep="","alt_samples_",alt)]]<-parList[[i]]$pars
alt95ResultsList[[paste(sep="","alt",alt)]]<-parList[[i]]$results
}
rm(parList)
save(alt95ResultsList,alt95ModList,file=paste(sep="","altResultsList_te95_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40
load("altResultsList_te95_20210805.RData")
sensVect<-numeric(length(altVect))
specVect<-numeric(length(altVect))
for(i in 1:length(altVect)){
sensVect[i]<-alt95ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-alt95ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
alt=altVect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-altVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nALT threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for ALT with diagnosis threshold te > 9.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on ALT (from 1 to 198 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
alt95ResultsList[[paste(sep="","alt",bestThr)]] %>%
kable(caption="Best (according to Youden's J) ALT model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 1.1684 | (0.3032,2.3067) |
| hazardous alcohol drinking (OR for specificity) | 1.0205 | (0.5890,1.4774) |
| suspected liver disease screening (OR for sensitivity) | 1.7509 | (0.7774,2.8940) |
| suspected liver disease screening (OR for specificity) | 0.1432 | (0.0992,0.1885) |
| hospital-based study (OR for sensitivity) | 1.0544 | (0.3193,1.9851) |
| hospital-based study (OR for specificity) | 0.7324 | (0.5285,0.9637) |
| random effects variance (logit sensitivity) | 0.3464 | (0.0635,0.8222) |
| random effects variance (logit specificity) | 0.9517 | (0.2286,2.0164) |
| overall sensitivity | 0.5919 | (0.4103,0.7667) |
| overall specificity | 0.8213 | (0.7763,0.8652) |
altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40
alt79ResultsList<-list()
alt79ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(altVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
alt<-altVect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_79==1,1,2)
datJags$y<-ifelse(dat$alt>=alt,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(altVect)){
alt<-altVect[i]
alt79ModList[[paste(sep="","alt_mod_",alt)]]<-parList[[i]]$model
alt79ModList[[paste(sep="","alt_samples_",alt)]]<-parList[[i]]$pars
alt79ResultsList[[paste(sep="","alt",alt)]]<-parList[[i]]$results
}
rm(parList)
save(alt79ResultsList,alt79ModList,file=paste(sep="","altResultsList_te79_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40
load("altResultsList_te79_20210805.RData")
sensVect<-numeric(length(altVect))
specVect<-numeric(length(altVect))
for(i in 1:length(altVect)){
sensVect[i]<-alt79ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-alt79ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
alt=altVect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-altVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nALT threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for ALT with diagnosis threshold te > 7.9kPa (F2).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on ALT (from 1 to 198 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
alt79ResultsList[[paste(sep="","alt",bestThr)]] %>%
kable(caption="Best (according to Youden's J) ALT model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 0.8833 | (0.3293,1.5962) |
| hazardous alcohol drinking (OR for specificity) | 0.8782 | (0.4385,1.4023) |
| suspected liver disease screening (OR for sensitivity) | 5.3853 | (2.6790,8.7592) |
| suspected liver disease screening (OR for specificity) | 0.0557 | (0.0352,0.0790) |
| hospital-based study (OR for sensitivity) | 0.9747 | (0.3977,1.6336) |
| hospital-based study (OR for specificity) | 0.6983 | (0.4467,0.9820) |
| random effects variance (logit sensitivity) | 0.8167 | (0.1439,1.8155) |
| random effects variance (logit specificity) | 1.2736 | (0.2933,2.8197) |
| overall sensitivity | 0.3709 | (0.2283,0.5160) |
| overall specificity | 0.9292 | (0.9019,0.9543) |
fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43
fib4115ResultsList<-list()
fib4115ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(fib4Vect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
fib4<-fib4Vect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_115==1,1,2)
datJags$y<-ifelse(dat$fib4>=fib4,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(fib4Vect)){
fib4<-fib4Vect[i]
fib4115ModList[[paste(sep="","fib4_mod_",fib4)]]<-parList[[i]]$model
fib4115ModList[[paste(sep="","fib4_samples_",fib4)]]<-parList[[i]]$pars
fib4115ResultsList[[paste(sep="","fib4",fib4)]]<-parList[[i]]$results
}
rm(parList)
save(fib4115ResultsList,fib4115ModList,file=paste(sep="","fib4ResultsList_te115_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43
load("fib4ResultsList_te115_20210805.RData")
sensVect<-numeric(length(fib4Vect))
specVect<-numeric(length(fib4Vect))
for(i in 1:length(fib4Vect)){
sensVect[i]<-fib4115ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-fib4115ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
fib4=fib4Vect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-fib4Vect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nFIB4 threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for FIB4 with diagnosis threshold te > 11.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on FIB4 (from 0 to 38.35 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
fib4115ResultsList[[paste(sep="","fib4",bestThr)]] %>%
kable(caption="Best (according to Youden's J) FIB4 model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 1.0804 | (0.1749, 2.3881) |
| hazardous alcohol drinking (OR for specificity) | 0.6838 | (0.3642, 1.0495) |
| suspected liver disease screening (OR for sensitivity) | 5.4163 | (0.9884,12.3368) |
| suspected liver disease screening (OR for specificity) | 0.2333 | (0.1094, 0.3800) |
| hospital-based study (OR for sensitivity) | 0.8995 | (0.1252, 2.0799) |
| hospital-based study (OR for specificity) | 2.2552 | (1.4912, 3.0942) |
| random effects variance (logit sensitivity) | 1.0580 | (0.1072, 2.7275) |
| random effects variance (logit specificity) | 0.6531 | (0.1213, 1.4914) |
| overall sensitivity | 0.5790 | (0.3201, 0.8425) |
| overall specificity | 0.8794 | (0.8378, 0.9191) |
fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43
fib495ResultsList<-list()
fib495ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(fib4Vect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
fib4<-fib4Vect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_95==1,1,2)
datJags$y<-ifelse(dat$fib4>=fib4,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(fib4Vect)){
fib4<-fib4Vect[i]
fib495ModList[[paste(sep="","fib4_mod_",fib4)]]<-parList[[i]]$model
fib495ModList[[paste(sep="","fib4_samples_",fib4)]]<-parList[[i]]$pars
fib495ResultsList[[paste(sep="","fib4",fib4)]]<-parList[[i]]$results
}
rm(parList)
save(fib495ResultsList,fib495ModList,file=paste(sep="","fib4ResultsList_te95_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43
load("fib4ResultsList_te95_20210805.RData")
sensVect<-numeric(length(fib4Vect))
specVect<-numeric(length(fib4Vect))
for(i in 1:length(fib4Vect)){
sensVect[i]<-fib495ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-fib495ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
fib4=fib4Vect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-fib4Vect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nFIB4 threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for FIB4 with diagnosis threshold te > 9.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on FIB4 (from 0 to 38.35 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
fib495ResultsList[[paste(sep="","fib4",bestThr)]] %>%
kable(caption="Best (according to Youden's J) FIB4 model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 0.7012 | (0.1610,1.4410) |
| hazardous alcohol drinking (OR for specificity) | 0.7515 | (0.4197,1.1372) |
| suspected liver disease screening (OR for sensitivity) | 4.7230 | (1.7261,8.3516) |
| suspected liver disease screening (OR for specificity) | 0.3127 | (0.1681,0.4883) |
| hospital-based study (OR for sensitivity) | 0.7404 | (0.1755,1.4538) |
| hospital-based study (OR for specificity) | 2.5483 | (1.7757,3.3484) |
| random effects variance (logit sensitivity) | 0.5324 | (0.0958,1.2093) |
| random effects variance (logit specificity) | 0.3937 | (0.0931,0.8414) |
| overall sensitivity | 0.6035 | (0.4076,0.7976) |
| overall specificity | 0.7950 | (0.7442,0.8447) |
fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43
fib479ResultsList<-list()
fib479ModList<-list()
cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)
parList<-foreach(i=1:length(fib4Vect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
fib4<-fib4Vect[i]
datJags<-list()
datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_79==1,1,2)
datJags$y<-ifelse(dat$fib4>=fib4,1,0)
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc
randSeed<-20210615
set.seed(randSeed)
jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
resultsModel4 <- resultsModel4 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}
parallel::stopCluster(cl)
for(i in 1:length(fib4Vect)){
fib4<-fib4Vect[i]
fib479ModList[[paste(sep="","fib4_mod_",fib4)]]<-parList[[i]]$model
fib479ModList[[paste(sep="","fib4_samples_",fib4)]]<-parList[[i]]$pars
fib479ResultsList[[paste(sep="","fib4",fib4)]]<-parList[[i]]$results
}
rm(parList)
save(fib479ResultsList,fib479ModList,file=paste(sep="","fib4ResultsList_te79_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43
load("fib4ResultsList_te79_20210805.RData")
sensVect<-numeric(length(fib4Vect))
specVect<-numeric(length(fib4Vect))
for(i in 1:length(fib4Vect)){
sensVect[i]<-fib479ResultsList[[i]]["overall sensitivity","Mean"]
specVect[i]<-fib479ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1
aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]
aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)
dfROC<-data.frame(
fib4=fib4Vect,
sensitivity=sensVect,
specificity=specVect,
YoudenJ=youdenJ
)
idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-fib4Vect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
geom_point(alpha=1,col="darkgrey") +
geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
geom_hline(yintercept=bestSens,lty=2,col="grey") +
geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nFIB4 threshold = ",bestThr)) +
labs(title=paste(sep="","ROC curve for FIB4 with diagnosis threshold te > 7.9kPa (F2).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on FIB4 (from 0 to 38.35 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
theme_light()
fib479ResultsList[[paste(sep="","fib4",bestThr)]] %>%
kable(caption="Best (according to Youden's J) FIB4 model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 0.4472 | (0.1478,0.8085) |
| hazardous alcohol drinking (OR for specificity) | 0.9615 | (0.5437,1.4109) |
| suspected liver disease screening (OR for sensitivity) | 5.3517 | (2.7271,8.4217) |
| suspected liver disease screening (OR for specificity) | 0.4250 | (0.2670,0.6050) |
| hospital-based study (OR for sensitivity) | 0.4106 | (0.1548,0.7384) |
| hospital-based study (OR for specificity) | 1.6224 | (1.1971,2.0830) |
| random effects variance (logit sensitivity) | 0.3803 | (0.0726,0.8606) |
| random effects variance (logit specificity) | 0.4677 | (0.1203,0.9764) |
| overall sensitivity | 0.7818 | (0.6616,0.8943) |
| overall specificity | 0.5663 | (0.5007,0.6301) |
tmp<-mod_te115_apri20 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
tmp %>%
kable(caption="Bivariate random effects model for APRI$\\geq2.0$ and te>11.5 (F4).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 1.7432 | (0.1091, 4.3214) |
| hazardous alcohol drinking (OR for specificity) | Inf | (0.1888, Inf) |
| suspected liver disease screening (OR for sensitivity) | 5.2212 | (0.5415,13.3010) |
| suspected liver disease screening (OR for specificity) | 0.2905 | (0.0032, 0.8721) |
| hospital-based study (OR for sensitivity) | 0.5133 | (0.0213, 1.4310) |
| hospital-based study (OR for specificity) | 7.9510 | (0.3947,24.2266) |
| random effects variance (logit sensitivity) | 1.2211 | (0.1405, 3.1288) |
| random effects variance (logit specificity) | 1.3605 | (0.0544, 4.5868) |
| overall sensitivity | 0.2876 | (0.0570, 0.5483) |
| overall specificity | 0.9918 | (0.9820, 0.9995) |
tmp<-mod_te95_apri20 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
tmp %>%
kable(caption="Bivariate random effects model for APRI$\\geq2.0$ and te>9.5 (F4).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 1.7008 | (0.1402, 4.0561) |
| hazardous alcohol drinking (OR for specificity) | Inf | (0.1927, Inf) |
| suspected liver disease screening (OR for sensitivity) | 9.0142 | (0.9858,22.4241) |
| suspected liver disease screening (OR for specificity) | 0.6857 | (0.0016, 2.1655) |
| hospital-based study (OR for sensitivity) | 0.5460 | (0.0182, 1.3568) |
| hospital-based study (OR for specificity) | 8.4439 | (0.2685,26.0366) |
| random effects variance (logit sensitivity) | 1.2412 | (0.1350, 3.1064) |
| random effects variance (logit specificity) | 0.7651 | (0.0649, 2.2119) |
| overall sensitivity | 0.1476 | (0.0270, 0.2852) |
| overall specificity | 0.9916 | (0.9818, 0.9990) |
tmp<-mod_te79_apri15 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
tmp %>%
kable(caption="Bivariate random effects model for APRI$\\geq1.5$ and te>7.9 (F2).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 0.9711 | (0.1073, 2.2673) |
| hazardous alcohol drinking (OR for specificity) | 10.8563 | (0.0749,21.1288) |
| suspected liver disease screening (OR for sensitivity) | 18.3817 | (2.7344,44.9702) |
| suspected liver disease screening (OR for specificity) | 0.1231 | (0.0050, 0.3097) |
| hospital-based study (OR for sensitivity) | 0.8613 | (0.1236, 1.8966) |
| hospital-based study (OR for specificity) | 5.9155 | (0.3831,16.7334) |
| random effects variance (logit sensitivity) | 1.3829 | (0.1927, 3.2764) |
| random effects variance (logit specificity) | 0.8537 | (0.0486, 2.5903) |
| overall sensitivity | 0.0732 | (0.0187, 0.1393) |
| overall specificity | 0.9916 | (0.9817, 0.9993) |
tmp<-mod_te115_gpr56 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
tmp %>%
kable(caption="Bivariate random effects model for GPR$\\geq0.56$ and te>11.5 (F4).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 2.0999 | (0.1663,5.1447) |
| hazardous alcohol drinking (OR for specificity) | 0.5082 | (0.1158,1.1127) |
| suspected liver disease screening (OR for sensitivity) | 3.1302 | (0.2450,7.9593) |
| suspected liver disease screening (OR for specificity) | 0.3137 | (0.0496,0.7137) |
| hospital-based study (OR for sensitivity) | 2.7554 | (0.1768,7.5154) |
| hospital-based study (OR for specificity) | 1.0458 | (0.3392,1.9581) |
| random effects variance (logit sensitivity) | 1.8676 | (0.1147,5.7481) |
| random effects variance (logit specificity) | 1.0584 | (0.0847,2.9765) |
| overall sensitivity | 0.2268 | (0.0333,0.4731) |
| overall specificity | 0.9845 | (0.9711,0.9956) |
tmp<-mod_te95_gpr56 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
tmp %>%
kable(caption="Bivariate random effects model for GPR$\\geq0.56$ and te>9.5 (F4).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 1.6621 | (0.1627, 3.8664) |
| hazardous alcohol drinking (OR for specificity) | 0.3858 | (0.0883, 0.8570) |
| suspected liver disease screening (OR for sensitivity) | 4.8646 | (0.8795,10.8084) |
| suspected liver disease screening (OR for specificity) | 0.3742 | (0.0331, 0.9331) |
| hospital-based study (OR for sensitivity) | 1.3581 | (0.2415, 3.0542) |
| hospital-based study (OR for specificity) | 0.8700 | (0.2400, 1.6434) |
| random effects variance (logit sensitivity) | 0.8703 | (0.0920, 2.3205) |
| random effects variance (logit specificity) | 0.8122 | (0.0805, 2.2122) |
| overall sensitivity | 0.2043 | (0.0534, 0.3795) |
| overall specificity | 0.9888 | (0.9788, 0.9969) |
tmp<-mod_te79_gpr32 %>%
as.data.frame() %>%
mutate(
Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")
tmp %>%
kable(caption="Bivariate random effects model for GPR$\\geq0.32$ and te>7.9 (F2).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
kable_styling(full_width = FALSE)
| Posterior mean | 95% HDI credible interval | |
|---|---|---|
| hazardous alcohol drinking (OR for sensitivity) | 1.6963 | (0.4509, 3.3217) |
| hazardous alcohol drinking (OR for specificity) | 0.2790 | (0.1215, 0.4660) |
| suspected liver disease screening (OR for sensitivity) | 6.1484 | (2.3796,10.9575) |
| suspected liver disease screening (OR for specificity) | 0.3898 | (0.1239, 0.7457) |
| hospital-based study (OR for sensitivity) | 2.4301 | (0.8859, 4.4199) |
| hospital-based study (OR for specificity) | 0.8137 | (0.4056, 1.2840) |
| random effects variance (logit sensitivity) | 0.6234 | (0.0888, 1.5131) |
| random effects variance (logit specificity) | 0.8072 | (0.0976, 2.0004) |
| overall sensitivity | 0.1888 | (0.0812, 0.3073) |
| overall specificity | 0.9676 | (0.9483, 0.9850) |